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A modification of an implicit approximate-factorization finite-difference algorithm applied to the two- 
dimensional Euler and Navier-Stokes equations in general curvilinear coordinates is presented for supersonic 
freestream flow about and through inlets. The modification transforms the coupled system of equations Into an 
uncoupled diagonal form which requires less computation work. For steady-state applications the resulting 
diagonal algorithm retains the stability and accuracy characteristics of the original algorithm. Solutions are 
given for inviscid and laminar flow about a two-dimensional wedge inlet configuration. Comparisons are made 
between computed results and exact theory. 


Nomenclature 



= constants defined in Eq. ( 8 a) 

A,B 

= matrices defined in Eq. (10b) 

A,B 

= Jacobian matrices defined in Eq. (10a) 

c 

= local speed of sound 

D 

diagonal matrix defined in Eq. (12b) 

e 

= total energy per unit volume defined in 
Eq. (2) 

E,F 

= vector arrays defined in Eq. ( 1 ) 

E,F 

= vector arrays defined in Eq. (5) 

f nf 2 Js 

= defined in Eq. ( 8 b) 

h 

= time step = Ar or Ar/2 

I 

= identity matrix 

f 

= Jacobian of transformation 

k^,ky 

= ox kyf(kl^-kl) 

mj,m2 

= defined in Eqs. (13) 


= defined in Eqs. (16b) 

M 

= viscous matrix defined in Eqs. (16a) 
and (16b) 

N 



= matrix defined in Eq. (13) 

P 

= dimensionless pressure, php^ 

Pr 

- Prandtl number 

Q*Q 

= vector arrays defined in Eq. (1) 

Re 

= Reynolds number 

R 

= right-hand side matrix defined in Eq. 

s; 

( 11 ) 

= viscous matrix defined in Eq. (15) 

/ 

= dimensionless time in physical space 

II 

-7 

II 

= matrices defined in Eq. (12c) 

U,V 

= Cartesian velocity component in x and 
y directions, respectively 

u,v 

= contravariant velocities defined in Eq. 
(4) 

x,y 

= physical coordinates 

yi*yi*yr^ 

= geometric derivatives 

a. 

= pl('l2c) 
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<X},d2,Ot^,Ot4 

= defined in Eqs. (16b) 

0 

= l/('/2cp) 

y 

= ratio of specific heat 


= central difference operators 

Ar 

= time step = /? 


= backward difference operators 


= implicit smoothing, usually 2 *f£ 


= explicit smoothing, usually = h 


= generalized coordinate normal to body 


= metrics defined in Eq, (3) 

6 

= defined as kj,u-^kyV 

K 

= coefficient of heat transfer 


= matrices defined in Eq, (12b) 

M 

= I/V 2 , or coefficient of viscosity 

f 

= generalized coordinate parallel to body 


= metrics defined in Eq. (3) 

P 

= dimensionless density, p/p^ 

T 

= dimensionless time in transformed 
space 

<t>^ 

= defined as 0.5 (7 - I)(w^ + ) 


= forward difference operators 


Introduction 

P roper design of the air inlets is a crucial factor in 
achieving the desired performance of supersonic 
airbreathing missiles and aircraft. In the typical inlet design 
case, mass flow adequate to the demands of the propulsion 
system with maximum total pressure recovery and minimum 
external drag is required for a range of freestream conditions 
and vehicle altitudes. Meeting these requirements has 
traditionally required extensive experimental development. 
Recent trends indicate, however, that certain features of the 
complicated inlet flowfield are amendable to computation, 
even though existing methods have some shortcomings. Rizzi 
and Schmidt* have applied a finite-volume approach to low- 
supersonic inlet flowfields with reasonable success. Arlinger^ 
and Reyhner^ have developed numerical methods that are 
based on the full-potential equation. Adoption of the full- 
potential equation for supersonic flow above a Mach number 
of 1.3 should be avoided, however, due to the difference in 
shock jump conditions used from the exact shock jump 
conditions; this procedure also disregards the total pressure 
loss in the inlet which is very critical in supersonic inlet design. 
Bansod^ and Hawkins et al.^ have applied explicit numerical 
methods to transonic flow about inlets. These methods, 
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however, even if extended to the supersonic regime, would be 
limited in convergence rate due to the necessary clustering o 
Doints near the cowl lip which reduces the integration step 
size More recently the first author has developed a method 
which is implicit in nature and faster than the previous 
methods. Other available methods’ * are for supersonic in- 


ternal flows. . i • 

This paper describes a new method of calculating the 
flowfield around and through two-dimensional inlets that 
does not have the limited Mach number range of potwtial 
methods; that is considerably faster than the methods of Refs. 
4-6- and that is valid for subcritical, critical, and supercritical 
mass flow rates. In the present method,’ a modification of an 
imolicit approximate-factorization finite-difference 
alRorithm'® " applied to partial differential equations is 
nresented which substantially reduces the total computation 
work. The modification takes the coupled system of equations 
into an uncoupled diagonal form which is easier to solve. The 
resulting diagonal algorithm retains the stability and many of 
the accuracy characteristics of the original algorithm. Since 
the governing equations written in generalized coordinates are 
cast in conservation law form, the complex features of the 
inlet flowfield are captured correctly. With the equations cast 
in generalized coordinates, mesh-generating routines are 
used to create the computational mesh. Clustering of grid 
points near solid surfaces is allowed. The implicit algorithm is 
then used to advance the unsteady equations in time with the 
ability of being able to take larger time steps than the previous 
explicit techniques, thus allowing faster convergence. 

Several different types of boundary conditions are 
required. At solid boundaries, the tangency condition is in- 
corporated as no flow through the boundary or a no-slip 
condition for the viscous cases. The upstream and lateral 
boundaries are specified so that freestream conditions may be 
specified (and these may be nonuniform). At outHow, it the 
flow is purely supersonic the variable are extrapolated. It the 
outflow is subsonic, as in a critical or subcritical calculation, 
the back pressure is specified and all other variables^ are 
calculated from characteristic-like equations of Kentzer. 


Governing Equations 

The partial differential equations governing the two- 
dimensional planar flow of an unsteady inviscid, nonheat- 
conducting, ideal gas can be written in nondimensional strong 
conservation law form” under the generalized independent 
variable transformation 

T=r, ^ = k(f.x.y). v = v(<.x,y) 


as follows: 


Equation (1) was integrated forward in time to a steady- 
state condition using an existing implicit Euler equation 
solver. * ^ 

The metrics in Eq. (1), ii, (x* are easily formed from 
the derivatives of , etc., using the relations 

riy^Jx^, ri, = -XrVx-yrVy ( 3 ) 

It is also convenient to define the velocities 

U=^,+i,u + i^v 


V=r,,+r,^U-i-r,yV W 

which are the so-called contravariant velocities along the £ 
and rj coordinates. Using these defined velocities, E and Fean 
be written in the compact form 


pU 


pV 

puU+kxP 

F^J-^ 

puV+rixP 

pvU-y^yP 


pvV+riyP 

{e+p)U-^,p _ 


{e+p) V-ri,p 


( 5 ) 


Note that once U and V are formed, the fiux vectors £ and F 
are not much more complex than E and F. To complete the 
problem, boundary and initial conditions must be specified. 


Boundary Conditions 

Along the body surface rj{x,yJ) =0 (the cowl, ramp, and 
inlet surfaces of Fig. 1), the condition of tangency in unsteady 
flow is enforced by 


F=owith 



( 6 ) 


For viscous now, t/=0 is used in Eq. (6) to produce no-slip 

conditions. u • j r 

The pressure on the body surface can be obtained trom ine 

normal momentum equation 


= iVxix'^^yVy)Pi + 0) 


q,+ [a,q + ijcE+^yF)/J]i + [v,Q + VxE + r,yF)/J],-0 ( 1 ) 

where 


" P 


pv 


pv 

pu 


pu^ +p 


puv 

E= 


F= 


pv 


puv 


pv^ +p 

e 


{e+p)u _ 


(e+p)v 


and the Jacobian 

In Eq. (1), p represents the pressure, p the density, u and i; 
the velocity components in the x and y directions, respec- 
tively, and e the total energy per unit volume. The following 
equation relates the pressure, density, and velocity com- 
ponents to the energy for an ideal gas 

e~p/{y — 7 ) -l-p -b ) /2 (2) 


where n is the direction normal to the body surface. 

At the internal outflow plane there are two possible How 
conditions. If the flow is supersonic, the conservative 
variables are zero-order extrapolated; whereas if the now is 
subsonic, the back pressure is specified and Kentzer s 
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Fig. I Surface geomelr>' of inlet. 
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approach is taken. In this approach, the characteristic 
compatibility relation, the projection of the momentum 
equations, and the energy equation are used to calculate the 
remaining independent variables. The boundary conditions 
for the subsonic outflow are 

yncw ^^Old =p^^<^ +Atfj 

(8a) 

P—Pb + By = const 

where 

/, = - (m + c ) +VxU,)-V{^yU(+rt^U^) 

_ /u + c\ 

-c(.kyV^ (^xPi +VxPv)-bv 

f2=-»iixV( +VxV^)-v{^yV^ +riyV,)-6/p 

fs= +VxP,,) +VyP,i) +bv (8b) 

and d and 6 are constants allowing for a linear distribution of 
back pressure across the outflow plane. 

At the external outflow boundary, the flow is supersonic 
and the conservative variables are extrapolated, as is usually 
the case. The inflow and outer boundaries are placed such that 
freestream conditions are specified. This completes the 
specification of the boundaries. 

Grid Generation 

The transformed equations are somewhat more com- 
plicated than the original Cartesian form, but offer several 
significant advantages. The main advantage is that boundary 
surfaces in the physical plane can be mapped onto rectangular 
surfaces in the transformed plane. Another significant aspect 
of the transformation is that grid points can be concentrated 
in regions that experience rapid change in the flowfield 
gradients. This is especially important in the present problem 
with numerous expansion and compression corners internal 
and external to the inlet. 

To take advantage of the generality of the transformed 
equations, one needs a fairly automatic method of generating 
a smoothly varying grid that conforms to arbitrary bodies and 
allows grid point clustering. The scheme that is chosen for the 
present application is the Thompson, Thames, and Mastin'^ 


method which has been altered by Sorenson and Steger’^ and 
further altered by the present authors. In this method the grid 
in the physical plane is defined by the solution of a Laplace or 
a Poisson equation. Grid points are arbitrarily specified on 
the body boundaries so that even if the Laplace equation is 
used, the generated grid is not orthogonal. The capability to 
select the location of boundary node points is one of the 
desirable features of the scheme and Eqs. (1) and (2) do not 
assume orthogonality. 


Numerical Algorithm 

An implicit numerical algorithm is used to solve the 
equations since in many flowfield problems it is desirable to 
take a larger time step than that permitted by a conventional 
explicit scheme. Such a situation may occur if the dependent 
variables experience a more rapid variation with space than 
with time. In addition, the unsteady form of the Euler 
equations [Eq, (1)] were solved to allow for regions of sub- 
sonic flow in the inlet to develop if a subcritical or critical 
flow condition existed. 

Conventional Form 

The basic numerical algorithm used was developed by Beam 
and Warming*^ and by Steger.^* It is second-order accurate 
in space and time, is noniterative, and is in a spatially factored 
form referred to as the “delta-form.'* A fourth-order dis- 
sipation term is appended to the right-hand side and in that 
location helps to control possible numerical instabilities. For 
either trapezoidal or Euler temporal implicit differencing, the 
delta form algorithm is given by 

+ (9) 

where for the convection terms and 5,, are second-order 
central difference operators, h — At or At/2 for first- or 
second-order two-level time differencing, and for convenience 
the spatial indices are deleted throughout. 

The Jacobian matrices AdinB are defined as 

B = (10a) 


where 


iy-S) , , ( 7 - 7 ) , 


>4 = 


— uv 


- (y-3)u 


“7w( J) + (7-7)w(w^ + i;^) t(^) ~ (iw^ + u^) 


0 0 

-(y-I)v (7-7) 

u 0 

-{y-l)uv yu 


B = 


0 

— uv 


iy-3) , , (7-7) , 


0 10 
V u 0 

-{y-l)u -(y-3)v iy-l) 

-(y-I)uv 7 ^^^ + (u^ +3v^) yv 


(10b) 
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Numerical dissipation terms are added to the implicit operators, with coefficient and the explicit side, with coefficient , to 
control nonlinear instabilities. Typically, e£ =0(A) and t/ =2f£. 

Diagonal Form 

A diagonal form of the implicit algorithm has been developed in Ref. 9. This algorithm retains many of the stabi ity and ac- 
curacy characteristics of the original scheme, but requires less computational work. A detailed derivation and analysis of the 
method can be found in Ref. 9. Here, we shall outline the basic development of the diagonal algorithm. 

The standard algorithm is rewritten in inviscid form for convenience as 

{/+ Wj/i " ) (/+ /i5, e" ) Ag = - A/ (5j + 5, F" ) = R " 

The similarity transformations which diagonalize /4 and B (see Warming et al. '*) are 

A = T^kji‘. b=t^a,t;' 

where 

U 0 0 0 

0 U 0 0 

0 0 U+cai+^y)''" 0 

0 0 0 f/-c(«i + £^)*^ 


kf=D[U.U.U+cai + ^^y)‘'^, U-cai + kj)'''’] = 




A, =£>[ K. K. v+c{ni+Vy) y-c(vi+Vy)'''‘] 

I 0 ot 

u kyP a(uAkyC) 

V -k„p a(v + kyC) 


Ty‘= -p-> (kyU-icyV) kyp^‘ 

(7-7)w] 

^(^^+c6) -i3[^^c+ (7 -/)u] 


anda=p/(V2c),/3 = l/(V2pc).e=^^« + ^^u, = [( 7 - 

l)/2] (u^ +u^),and, for example, ky=kyl (k^ + kl) 

Relations exist between and of the form 

= (i3) 

where 

~] 0 0 O' 

0 rrij jfini 2 

0 jiim2 + 

whh rrjf = (^xVx'^iyVy)* (^xVy‘~‘iyVx)* ft=l>/2, 

where y = 1 for the matrix N and y= - 1 for the inverse matrix 

N-'. 

Applying the similarity transformations [Eqs. (12) and (13)] 
to Eq. (II) and factoring and out of the spatial 
operators, we have the diagonal algorithm 

r^(/+A5^ApN(/+/i6,A;)(r;M"A4=^" 04) 


The new implicit operators, (7-l-/t6{A|) and (/+/?5,,A,j) , 
are still block tridiagonal operators but now the blocks are 
diagonal in form such that the equations can be reordered into 
four separate scalar tridiagonal operators. This has a large 
positive impact on the solution process discussed below. 

The solution process for the implicit part of Eq. (14) 
consists of: 1) S/ = ( Tf M a matrix-vector multiply at 
each grid point, since is known analytically; 2) four 
scalar tridiagonal inversions for the operator 
S, = [ /+ A^ ] ' Sj: 3)S, = N - ' S, , ^ a matrix-vector 

multiply at each point; 4) S^= [/-l-/i6^A, [ S^, four more 

scalar tridiagonal inversions; 5) A^=F"S,^, another set of 
matrix-vector multiples; and finally 6) to 

update the solution. This contrasts with the two block 
tridiagonal inversions required in Eq. (1 1). 

An operation count for the diagonal form of the implicit 
algorithm yields 233 multiplies, 125 adds, and 26 divides, 
totaling 384 operations. For the standard algorithm in the 
transformed coordinates, the operation counts are 410 
multiplies, 326 adds, and 10 divides for a total of 746 
operations. The use of the diagonal algorithm produces a 
33% savings in computer time on a CDC 7600 for a realistic 
calculation. 
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Viscous Form 

The viscous form of the equations and algorithm has been 
described in detail by Steger** along with a discussion and 
justification of the “thin-layer” model. In the thin-layer 
approximation, viscous terms along the body, f direction, are 
neglected and those in the or near normal direction, are 
retained. 

The “thin-layer” Navier-Stokes equations in general 
curvilinear coordinates are 


(15) 

where 


0 


3,s=a,7- 


li{ril+riy)v^ + (ti/ 3 )rjy{rj^u^ +riyV^) 

[kPr-'{y-]) 

+ li{vi+Vy) (u^+v^)^/2 

+ iti/3){ri^U + riyV){Vj,U^+r]^V^)] 


The viscous form of the conventional implicit algorithm is 
(I+h6^A" V ^A^J) (I+hSy^B" '7,A,J 

-hRe-'b„M’')Aq’' = -At(bf^E’' +5,F'’ -Re-'bJ”) 
-efy-'[(V{Aj)^ + (V,A,)^]7<7'' (I6a) 

where 


0 




ni2j 




0 0 0 ~ 

oc,d,(p-') a;3,(p-') 0 

a;3,(p"') aj3,(p‘') 0 

f^42 ^43 ^44 _ 


m2i = -otid^(u/p) -a2d^(v/p) 
mjj = -ot2drf{u/p) -asd^iv/p) 

m4i [ - {e/p^) -h (w^ + y^)/p] 

-ajd^ (u^/p) - 2 a 2 d^ (uv/p) /p) 

ni42= -ot^d^ (M/p) -nj2j 


^43 = 

ni44 = oi4^^(P'^) 

cxj =p[ ct2 = 

Gi3 =pWx-^(^/3)nl]> ot4='yKPr’^ ( 1 ^^) 

Here Re is the Reynolds number, Pr the Prandtl number, p 
the dynamic viscosity, and k is the coefficient of thermal 
conductivity. 

The inclusion of the viscous term in the second implicit 
operator of Eq. (16a) makes it difficult to apply the 
diagonaljzation to the viscous scheme. This is because the 
matrix M does not have the same eigenvector matrices as B 
and therefore cannot be simultaneously diagonalized. 
Tannehill et al. suggest a procedure which will allow us to 
use the diagonal form. They suggest that the viscous terms in 
the implicit operators can be neglected for steady-state 


problems without affecting the stability or accuracy of 
moderate to high Reynolds number flows. It has been our 
experience that this can be done as long as a nonzero value of 
implict dissipation coefficient €/ is used. The elimination of 
the viscous terms in the implicit operators along with the 
diagonal algorithm produces a significant reduction in 
computational cost for steady-state problems. 

Results 

The results presented here demonstrate that the present 
numerical technique can calculate complicated two- 
dimensional inlet flowfields for varying conditions. Also 
discussed are the advantages of using the diagonal algorithm 
and of neglecting the implicit viscous term over the use of the 
standard algorithm. In all cases, the grid is initially calculated 
by a separate program and is kept nontime varying 
throughout the calculation. 

A two-dimensional cowl-ramp inlet system with a 10 deg 
ramp and 20 deg wedge-cowl is used as the test case at a 
freestream Mach number, A/o^=2.0. Different internal 
outflow conditions are used to produce supercritical and 
subcritical flows for inviscid conditions. 

A laminar viscous supercritical case is also presented. Since 
there were no experimental data available, the solution ob- 
tained was not optimal with respect to the grid point 
distribution. The intention here is to demonstrate the ability 
of the present algorithm to solve viscous as well as inviscid 
flows. A typical grid system is presented in Fig. 2, where grid 
lines can be clustered next to the body surfaces either to 



Fig. 2 Grid created by mesh generation program. 
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Mx = 2.0 ■'Theory “ ° 

‘^DcoWL " ° “numerical = 0.81233 


NUMERICAL RESULT 

2-D SHOCK THEORY 



Q I 1 1 t 1 J I I I I 1— —I 1 1 * 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 


SURFACE OF RAMP - DIFFUSOR - CHANNEL 



Fig. 4 Comparison between inviscid numerical results and two- 
dimensional wedge theor>'. 

M, = 2.0 PbaCK = 8.2 MtHEORY = °-814 
Cd^Owl " 0.178221 M^uMERICAL = 0.820 


SYMBOL 



Fig. 5 Pressure variations of subcritical inviscid results showing 
normal shock wave at the throat. 

improve the accuracy of boundary conditions or for viscous 
resolution. 

The numerically generated flowfield as viewed by density 
contours for a supercritical case is presented in Fig, 3, In the 
supercritical flow situation, there is supersonic flow 
throughout the flowfield, causing numerous shock waves and 
expansion waves which form a complex interaction system 
within the inlet channel. Figure 4 shows a comparison be- 
tween the numerical results and the theory for the pressure 
along the solid surfaces. Since the two-dimensional shock 
theory does not account for interactions of the waves there are 
slight differences in the shock locations between the 
theoretical and numerical solutions. 

A subcritical case is presented in Fig. 5 w'here, based on 
theoretical calculations, a back pressure of plp^ = S.2 is 
imposed. Shown are the pressures along the upper and lower 
surfaces and also along the middle coordinate line of the inlet. 
A strong shock wave located at the throat is considered to be a 
product of a critical to subcritical case in near-design con- 
ditions. This shock w^ave developed from a compression wave 
which formed at the internal outflow plane and then moved 
upstream. The normal shock wave sits just outside the dif- 
fuser throat causing a slight spillage of flow' onto the external 
cow l surface. This increases the cowl drag. 


= 2.0 “theory = ®-814 

CojjowL ■ ® “numerical = °-808 


NUMERICAL LAMINAR 

2-D SHOCK THEORY 



SURFACE OF RAMP - DIFFUSOR - CHANNEL 



0 I I I I I -i 1 1 1 1 1 1 ' ' ' 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

INTERNAL SURFACE OF COWL 
Fig, 6 Comparison between numerical laminar results and two 
dimensional wedge theory’. 



Fig. 7 Pressure contours for supercritical viscous laminar flow, 

A/oo=2.0. 


A viscous laminar calculation for the supercritical flow 
conditions at a Re/h of 10^ is shown in Figs. 6 and 7. In Fig. 
6 the pressure along the inlet surfaces are compared with 
theory. Even though viscosity dissipates the effects of the 
shock wave, the overall solution is comparable to the inviscid 
case. A small reverse flow bubble occurs on the expansion side 
of the internal wall of the ramp. Pressure contours in the cowl 
tip region (Fig. 7) show the shock wave that was formed. 

The supercritical case shown in Figs. 3, 4, 6, and 7 is used to 
investigate the effect of the diagonal algorithm, the dropping 
of the implicit viscous term on the convergence history, and 
the accuracy and efficiency of the numerical calculations. The 
conventional algorithm is used as the reference case. All 
comparisons are for 900 iterations (which assures con- 
vergence) and were performed on a CDC 76(X) computer for a 
66 X 36 grid point mesh system. 
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Fig. 8 Convergence history of diagonal 
algorithm as compared to standard algorithm. 



The inviscid calculation using the conventional algorithm 
requires 794 s of CPU time, whereas the diagonal algorithm 
requires 522 s, a 34^o savings. The convergence history for the 
two calculations are shown in Fig. 8. The residual is the root 
mean square of the right-hand side of Eq, (9) and can be 
considered as the norm of R'^ . Both cases converge at the 
same rate and reach identical steady-state solutions. 

The viscous calculation requires 975 s of CPU time for the 
conventional algorithm. By eliminating the viscous term in the 
implicit operator, the time is reduced to 834 s, a 14.5*70 
reduction. The converged solution is identical to the first 
solution. Furthermore, combining the viscous approximation 
with the diagonal algorithm brings the run time down to 567 s, 
or a 42*Vo reduction over the conventional scheme. In all three 
cases the convergence history and converged solutions are 
almost identical. 
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